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ABSTRACT 

In Paper I, we presented a one-dimensional hydrodynamic model for the winds 
of close-in exoplancts. However, close-in cxoplancts are tidally locked and irra- 
diated only on the day sides by their host stars. This requires two-dimensional 
hydrodynamic models with self-consistent radiative transfer calculations. In this 
paper, for the tidal locking (two-dimensional radiative transfer) and non-tidal 
locking cases (one-dimensional radiative transfer), we constructed a multi-fluid 
two-dimensional hydrodynamic model with detailed radiative transfer to depict 
the escape of particles. We found that the tidal forces (the sum of tidal gravity 
of the star and centrifugal force due to the planetary rotation) supply significant 
accelerations and result in anisotropic winds. An important effect of the tidal 
forces is that it severely depresses the outflow of particles near the polar regions 
where the density and the radial velocity are a factor of few-ten smaller than 
those of the low-latitude regions. As a consequence, most particles escape the 
surface of the planet from the regions of low-latitude. Comparing the tidal and 
non-tidal locking cases, we found that their optical depths are very different so 
that the flows also emerge with a different pattern. In the case of the non-tidal 
locking, the radial velocities at the base of the wind are higher than the merid- 
ional velocities. However, in the case of tidal locking, the meridional velocities 
dominate the flow at the base of the wind, and they can transfer effectively mass 
and energy from the day sides to the night sides. Further, wc also found that the 
differences of the winds show middle extent at large radii. It expresses that the 
structure of the wind at the base can be changed by the two-dimensional radiative 
transfer due to large optical depths, but the extent is reduced with an increase in 
radius. Because the escape is depressed in the polar regions, the mass loss rate 
predicted by the non-tidal locking model, in the order of magnitude of 10^° g 
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s"^, is a factor of 2 lower than that predicted by one-dimensional hydrodynamic 
model. The results of tidal locking show that the mass loss rate is decreased 
a order of magnitude, only 4.3 x 10^ g s~^, due to large optical depths on the 
night side. We also found that the distributions of hydrogen atoms show clear 
variations from the day side to night side, thus the origin of the excess absorption 
in Ly a should be reexamined by using multi-dimensional hydrodynamic models. 

Subject headings: hydrodynamics - planetary systems - planets and satellites: 
atmospheres - planets and satellites: individual (HD 209458b) 
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1. INTRODUCTION 

The excess absorption in Lyman a. of HD 209458b and HD 189733b have been 
found by Vidal-Madjar et al. (2003) and Lecaveher des Etangs et al. (2010). Prom the 
excess absorption those authors suggested that the hydrogen clouds of HD 209458b and 
HD 189733b are above their Roche lobe. Most of the explanations attribute the excess 
absorption to mass loss from the surfaces of exoplanets (Lammer et al. 2003; Lecaveher 
des Etangs et al. 2004; Yelle 2004, 2006; Tian et al. 2005; Garcia Muhoz 2007; Penz et al. 
2008; Murray-clay et al. 2009; Lammer et al. 2009; Guo 2011 (Paper I)). To simphfy the 
calculations, all the hydrodynamic models above are one dimensional. Recently, two- and 
three-dimensional hydrodynamic models have been presented by Stone & Proga (2009) and 
Schneiter et al. (2007), however, some challenging physical processes have been neglected 
by them. All the above one-dimensional models described thermal particle escape, namely, 
the particles of the planets escape from their gravity potential wells as a result of the 
irradiation from the host star. However, the charge exchange between the stellar wind 
and the planetary escaping exosphere can result in the same observation phenomenons 
(Holmstrom et al. 2008). The loss of nonthermal neutral atoms via interactions between 
the stellar wind and the exosphere have been discussed by many researchers (Erkaev et al. 
2005; Holmstrom et al. 2008; Ekenback et al. 2010). Their results showed that the fresh 
neutral atoms produced by the charge exchange can match the Ly a excess absorption. 

Although current theoretic models predicted the mass loss at the upper atmosphere 
of Hot Jupiters, no clear observational evidences support the thermal or non-thermal 
explanations. Either the energetic HI of stellar origin or thermal HI populations in the 
planetary atmosphere can fit the observations of Ly a (Bcn-Jaffel & Hosseini 2010). 
Moreover, Koskinen et al. (2010) found that the excess absorption can be explained solely 
by absorption in the upper atmosphere and the process of charge exchange might not be 
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necessary. According to the examples from our solar system, however, the charge exchange 
can play an important role in explaining the excess Ly a absorption (see Lammer 2011 for 
more details). 

The extent of understanding the excess absorption depends strongly on the physics 
described in theoretical models. To fully distinguish the thermal or non-thermal processes, 
the single-fluid and one-dimensional models are not enough. Yelle (2004) and Garcia 
Munoz (2007) have used the diffusion approximation to depict the velocity differences 
between different species, and photochemistry, ionization and recombination have also been 
included. Moreover, Guo (2011) presented a multi-fluid one-dimensional model to depict 
the atmospheric escape and predicted reasonable mass loss rates for HD 209458b and HD 
189733b. The multi-fluid model is more advantageous than single-fluid one because the 
interactions among different species between the stellar and planetary winds are depicted 
by its own continuity, momentum and energy equations. However, the model of Paper I can 
be improved in two points at least. First, the one-dimensional hydrodynamic model only 
simulates the case along the line connecting the center of planet and the center of the star. 
For the close-in expolanets, the tidal forces between the star and planet tend to circularize 
the orbit and synchronize the rotation of the planet with the orbital period (Guillot et al. 
1996; Trilhng 2000). Thus, for the close-in exoplanets it is reasonable to assume that their 
same sides are always facing the star. If the irradiation and tidal effects of the host star 
are included, multi-dimensional models are more suitable. Second, the constraint of stellar 
wind has been omitted. Stone & Proga (2009) discussed the interaction between the stellar 
and planetary winds and found that the planetary wind is confined to a small volume by 
the stellar wind. However, their model used the assumption of a single-fluid so the processes 
of charge exchange could not be simulated. 

These issues as above motivated us to construct a two-dimensional hydrodynamic 
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model with realistic treatment of radiative transfer and tidal forces. We have developed a 
time-independent model with the one-dimensional assumption, however, it is impractical 
to extend the method to two-dimension because of the existence of critical points. In this 
paper we constructed a time-dependent two-dimensional model including the micro-physics 
processes, the influence of tide and radiative transfer. The interactions between the stellar 
and planetary winds will be discussed in the further work. 

This paper aims to calculate the escape of atomic hydrogen and protons through the 
solution to their mass, momentum and energy equations. The microscopic physics processes 
are covered in the mass, but not in the momentum and energy equations (Section 2.1). The 
processes of radiative transfer are discussed in Section 2.2. We used semi-implicit scheme to 
calculate these equations (Section 2.3) and tested the model in Section 3. We incorporated 
one- and two-dimensional radiative transfer into the hydrodynamic models (Section 4), 
and focus on the differences between them. In Section 5 we discussed the influence of the 
boundary conditions and stellar activity. Our results are summarized in Section. 6 

2. THE MODEL 

2.1. Equations and Assumptions 

In Paper I (Guo 2011) we described a steady-state, radial expansion of plasma 
containing three species: atomic hydrogen (h), protons (p) and electrons (e). In this paper, 
we extend the previous work to two-dimensions. But thermal electrons are not included. 
Atomic hydrogen (h) and protons (p) have their own continuity equations and are described 
by a particle density rih and Up. Since this model deals with a mixture of atomic hydrogen 
and protons, the following processes are considered: photoionization and recombination. 
Further, we assumed that the velocities and temperatures of the particles are the same. 
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namely, Uh — Up — u and Th — Tp — T. The assumption is acceptable because the mass loss 
rates predicted by this model are higher than 10^ g (see Section 3 and 4), if the mass 
loss rates are below the value the velocities and temperature should be calculated separately 
because of the decoupling of particles (Guo 2011). As in Paper I, we did not include H2 in 
this model because the thermosphere of close-in planets should be composed primarily of H 
and H+. The location of transition from H2 to H is about l.lRp (Yelle 2004). 

We described the planetary wind by using multi-particle two-dimensional hydrodynamic 
equations that can be written as 

+ V • (n,u) = 5, (1) 



dt 



dinn) , , 

+ V • (nuu) + Vp = n^^t (2) 
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^ +V-[nu(e + y)] + V(pu) = naea:t-u + g, (3) 

where nj{j—h,p) is particle number density of hydrogen atoms and ions. n—Uh+np is the 
number density of the g whole. Sj denotes the production/destruction of particles. 

The specific internal energy is e = and = + Uq is the sum of squares of velocities, 
where f and 6 velocities are u^- and uq. 7 is set to 5/3. The accelerations produced by 
the gravities of planet and star as well as the centrifugal forces due to the rotation of the 
planet around the star are included in the term of ai.ext- The net heating rate is expressed 
as Q — H — L, where H is the heating from irradiation of the star and L is cooling due to 
the atmospheric radiation. 

The hydrodynamic equations are solved in spherical polar (r,^) coordinates. Here 
^=0 is the direction towards the star, and O—tt is the direction away from the star. In 
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conservation form, these equations can be expressed as 

dt dr do 



(4) 



where Q is the state vector, F and G are the flux vectors, 5*1 denotes the contributions 
from the planetary gravity and the tidal forces of stars, and S2 includes the source/sink 
terms that are relative with photoionization/recombination and heating/cooling. 

The state vector Q are defined as 
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The flux vector [F, G] and the source terms 5*1 are deflned as 
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where a = \/TZT is the isothermal sound velocity, and and are the accelerations in f 
and ^ directions, respectively. The accelerations are given by 



Or 



III' 
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For the position vector x = (r, 9), the effective potential U is given by 
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C/(x) 



(11) 



|x| Ix-Dl 2' 

where D is the distance of star, fl = [G(M* + Mp) / D^Y^'^ is the angular velocity of orbital 
plane, and and Mp are the stellar and planetary mass, respectively. Finally, we can 
obtain the accelerations as 



GM„ , ^^^.(Z^cos^^-r) ^.^ ^. 

+ TT^TT- — ^ — ttt: — G— COS 9{D — r cos 9) 



(i^2 + ^2 _ 2L)r cos ^)3/2 L)3 



(12) 



GM^Dsin^ 



+ G'— ^sin^(L)-rcos^). 



(13) 



(/)2+^2_2DrCOSe)3/2 D3 

In Equation (12), the first term represents the gravity of planet. The second and third 
terms as well as the two terms of Equation (13) are the accelerations produced by the 
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gravity of host stars and orbital motion, respectively. Except the first term of Equation 
(12), we refer to the other terms as "tidal accelerations". How do the tidal accelerations 
alter the behaviors of flows? We can estimate this by considering some special cases. In 
the case of 9 =0°, Equation (12) can be simphfied to = + - G^{D - r). 

This form is consisitent with the Equation (15) of Garcia Munoz (2007) and hints that 
the tidal forces supply positive acceleration along the line connecting the planet to the 
star. However, the radial tidal acceleration gradually decreases with the increase of 9 and 
become negative at 60° ^ < 120°. This means that the escape of particles in the polar 
regions can be depressed by the tidal forces. In addition, we also note that dg is always 
negative (positive) if the 9 is smaller (greater) than ~ 90°. These characteristics will also 
be discussed further in Section 4. 

The Coriolis forces have been neglected in the calculations. In the regions close to 
host stars, the escape velocities of particles which are in the order of magnitude of 20-40 
km s~^ are smaller than the orbital velocities (in the order of magnitude of 100 km s~^). 
This means that the planetary wind can be deflected by Coriolis forces. However, our 
calculations are limited to the regions close to exoplanets where the velocity of the wind is 
much less than the orbital velocity. Thus, Coriolis forces can be neglected safely near the 
planet. 

The source terms 5*2 can be expressed as 

r ^ sin 6 T^^c 

Qllfpho J.2 sine Tree 

-^2= , (14) 



sin9{H — L)/mH 

where m.H is the mass of hydrogen atom, 7pho = Ylv ^"hJ" (where F^, is photon energy 
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ed, frequency u, Ti, is the optical depth at frequency v, the cross section of photoionization 
of hydrogen is (t^ = 6 x lQ-^^{hv /UMV)-^ cm^ (Murray-Clay ct al. 2009)) and 
7rec = 2.7 X 10~^^( j^)~°-^ cm~^s~^ (Murray-Clay et al. 2009) represent the photoionization 
and recombination rates, respectively. In the continuity equation, the sources and sinks to 
the particle flux density are due to photoionization and recombination. In energy equation, 
the sources and sinks are heating and cooling, respectively. 

The heating processes in the mixed flows are complex. The models of Murray-Clay et 
al. (2009) and Paper 1 in which the heating from the host stars is simplified to a single 
photon of 20 ev predicted reasonable mass loss rates. However, the assumption of single 
photon energy could miss some important physical processes, such as the penetration of 
higher energy photons deeper in the atmosphere. In fact, the stellar radiation ionizes the 
species to produce high-energy photoelectrons that share their energy with other species 
via collisions. The heating of exoplanet atmospheres mainly comes from photoelectrons 
produced by photoionization, therefore, here we included EUV irradiance shortward of 912 
A, namely, the threshold value of photoionization in hydrogen atoms. 

In an ionized atmosphere composed of H and H+, photoelectrons distribute their 
energy via the collision with thermal electrons. H+ obtain energy from thermal electrons via 
Coulomb collisions while hydrogen atoms share the energy of thermal electrons by elastic 
and inelastic collisions. Here we assumed that H and H+ have the same temperatures, 
therefore, the total heating of the atmosphere can be written as 



For different photon energies, we defined the net heating efficiency as r^,^ = {hv—l?>.Q)/hv 
that is the maximum heating efficiency. In fact, the fraction of photon energy deposited as 
heat in the atmosphere of exoplanet should be smaller than the definition above. It also 
shows that the heating efficiency is higher for the photons with higher energy. Equation 




(15) 
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(15) hints that the contributions of heating by lower energy photons are comparable with 
those of higher energy photons because the heating of atmosphere is proportional to heating 
efficiency, but is inversely with cube of cross section of photoionization. 

The cooling rate adopted in the models of Murray-Clay et al. (2009) and Paper I 
assumed that the atmosphere is optically thin in the Ly a line. Murray-Clay et al. (2009) 
have also validated that a photon can diffuse out of the wind in a typical mass loss rate 
of 10^° g s~^ (see Appendix C of Murray-Clay et al. 2009). However, the observations 
indicate that the atmosphere is optically thick even in the wings of of Ly a. Therefore, it is 
uncertain if the Ly a. cooling should be included. Here we included the radiative cooling by 
recombination and free-free radiation, namely, 

L — Lj-ec + Lffeicgcm''^s'~^. (16) 

In addition, we can not calculate the infrared emission of Hj as done by Yelle (2004) 
and Garcia Munoz (2007) due to the lack of photochemistry in our model. The cooling of 
H3 could be important because it can modify the location of maximum temperature, but 
the mass loss rates are not affected by the cooling of Hg (Penz et al. 2008). 

There are not resolved emission spectra from the host stars of exoplanets. For HD 
209458, we used the EUV spectra of the sun as its emission spectra. Those data are taken 
from Richards et al. (1994). Our calculations are completed with a solar proxy P10.7 = 80 
that corresponds to a middle-low level of stellar activity. 

2.2. Radiative transfer 

The radiative transfers are the most important processes in determining the atmospheric 
structure of planet. In the case of non-tidal locking, both hemispheres of the planet can 
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receive the irradiation from its host star. Thus, the one-dimensional radiative transfer can 
approximately depict the process. In the assumption of ID radiative transfer, the optical 
depth at frequency v can be written as 



Clearly, Equation (17) shows that the optical depth is function of radii, but is not 
relative with Q. Thus, it cannot correctly depict the penetration of the EUV radiation from 
the face of day-side to night-side if the planets are tidally locked by their host stars. To 
calculate an accurate optical depth in the case of tidal locking, two-dimensional radiative 
transfer calculations are needed. Here we present a simplified treatment in which both 
the impact factor q{r,9) = rsin^ and latitude 9 determine the final optical depths. For 
each grid point, the nodal points at every radii are found along the opposite direction of 
stellar radiation, and the values of latitude at those radii are determined (see Figure 1). 
For a given grid point P(i, j), the nodal point at the next radius rj+i can be expressed as 
P(i-|-1, j') (Note that the impact factor q{r,6) of the given grid point P(i, j) is same with 
that of P(i-I-1, j')-)- Generally, the latitude of P(i-I-1, j') is between grid points 9j and Oj-i 
rather than consistent with the grid points of angular direction. Finally, the corresponding 
number density of P(i-I-1, j') at radius rj+i can be approximated as fih = +"'+i>j-i . 
Therefore, the optical depth along the z-direction (the direction of stellar radiation) can be 
approximated to 



The optical depth is set to 10^° when the impact factor q is smaller than Rp. Equation (18) 
is used to calculate the radiative transfer of tidal locking. 

The simplified calculations represent the variations of the optical depth at all grid 




(17) 




(18) 
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points. We believe that small differences due to the approximation do not affect the results 
quantitatively. 

2.3. Numerical Method and Boundary Conditions 

The two-dimensional time-dependent multi-fluid model is integrated in time by using 
the fourth-stage Runge-Kutta method. For the space derivative, central difference scheme 
is used. To solve Equation (4), we must handle those source terms in the continuity and 
energy equations because of their stiffness. An explicit method could result in a long 
time-marching and incorrect results. Thus, we choice a semi-implicit method to deal with 
S2 of Equation (4). 

Equation (4) is discretized on a grid of 221 radial and 121 angular cells. We selected 
exponentially spaced grid in radii, 

n^gn_i,g^exp( ^^^^^_^ ) (19) 

where i is a number of a given grid point, IMAX is the maximum number of grid points, 
and ri and tjmax denote the first and last grid points. The grid size increases exponentially 
with the increase of radius. We used a uniform grid in 9 direction. 

The lower boundary conditions at r — Rp are fixed {Rp is the planetary radius). The 
effective temperatures of hot Jupiter estimated by Burrows et al. (2001) are in the range of 
1000 - 1500K, which is consistent with the IR thermal emission measurements (Deming et 
al. 2005; Charbonneau et al. 2005; Deming et al. 2006). Thus, we set T = 1500K at the 
bottom of the flow. Here the pressure at the base of the wind is maintained on the level of 

= 1 dyn cm^^. The number density of H+ at the lower boundary is set to Hp = lO^'^nh- 
The calculation results are insensitive to the ratio of Up/uh- Finally, we assumed uq — Q at 
the bottom of the atmosphere. 
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We selected tjmax = 7Rp as the upper boundary, which approximately denotes the 
stand-off distance of the planetary wind due to the interaction with the stellar wind. At 
the upper boundary, we used the outflow boundary conditions. For supersonic outffow, the 
boundary condition is exact (Stone & Norman 1992). In this paper, we only consider the 
case of supersonic outflow. The planetary winds can be depressed to subsonic if the stellar 
winds are strong enough (Garcia Munoz 2007; Murray-Clay et al. 2009). A full description 
to the process is beyond the scope of this paper. We will discuss the effect of the interaction 
with the stellar wind in further work. 

We are interested in steady-state solutions. For all hydrodynamic solutions, the 
steady-state is reached when the relative change in the conservative variables from one time 
level to the next drops below 10^*^. We use a normalized measure defined by (Toth et al. 
1998) 



A2Q 



where iVvar is the number of conserved variables Qi, and the superscripts indicate the 
number of time level. 



3. TEST THE ONE-DIMENSIONAL MODEL 

To test our model, we started from ID planetary winds and applied the ID model to 
a typical and particular planet sample: HD 209458b. We take the the radius Rp = lARj 
and the mass of Mp = 0.7Mj, where Mj and Rj are the mass and radius of the Jupiter. 
The semi-major axis a is set to 0.047AU. The mass of host star is M^, = 1.148M0 



(http://exoplanet.eu/), where Mq is the mass of the sun. Note that the tidal acceleration 



in the ID model is written as 



iD^-^D^^^-'^^ (21) 
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which is same as that of Garcia Munoz (2007). 

Yelle (2004), Tian et al. (2005), Garcia Mufioz (2007) and Penz (2008) have developed 
time-dependent one-dimensional hydrodynamic models. Thus, a direct comparison with 
their results can vahdate the rehabihty of our model. 

With the stellar and planetary physical parameters, the time-dependent results are 
shown in Figure 2. It is clear from Figure 2 that the time-dependent model predicts the 
same trends of temperature, velocity and particle number density as do the other models. 
Those studies that the temperature can attain 8000-lOOOOK at 1.5ii!p are in agreement with 
our results. The velocities predicted by our model are higher than those of Yelle (2004), 
but are similar with those of Tian et al. (2005), Garcia Mufioz (2007) and Penz (2008). 
The fact that as much as 50% of hydrogen at r = 2.Rp is ionized is a consequence of being 
irradiated by the host star. Due to lack of photochemistry, the ionization state cannot be 
predicted accurately by our model. In fact, the atmosphere can be thought of as if made 
up solely of hydrogen if the amounts of heavy constituents are small enough. 

The total mass loss rate of 3.2 x 10^° g predicted by our model is comparable with 
the those of Yelle(2004), Tian et al. (2005), Penz (2008) and Garcia Mufioz (2007) (In 
Table 1, we summarized the cases with dyn cm~^ and P10.7 = 200. The mass loss 

rate of 2.2 x 10^^ g s~^ is comparable with that of Garcia Munoz (2007)). Some chemical 
processes are not included in our model. Thus, the net chemical expense of energy cannot 
be calculated. However, the escape of particles in the atmosphere of planet should be driven 
mainly by irradiation of host star. Our results show also that the mass loss rates are in 
good agreement with those of others, thus the effect of chemical energy could be slight. 

Yelle(2004) and Garcia Muiioz (2007) have found that the escape rate is energy limited, 
namely, the escape rate varies roughly with the stellar flux. To test the effect of stellar 
flux, we multiplied the solar spectrum by factors of 2.2 and recalculated the mass loss rate. 
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We found that the mass loss rate predicted by the model of 2.2 times solar spectrum is a 
factor of 3.7 higher than that of solar spectrum. Similar conclusions have also been found 
by Penz et al. (2008). From Table. 1 in Penz et al. (2008), it is clear that the ratio of mass 
loss rates are not proportional to the stellar flux. The difference can be explained by a few- 
reasons. First, photochemistry has been neglected in our model and Penz et al. (2008). 
However, it is unclear whether the photochemistry can affect the ratio of mass loss rates. 
We will discuss this in future work. Second, we note that the net heating efficiency in our 
model is defined ets rju — {hi/ — 13.6) /hv, which varies with the frequency, but not with 
radii. In fact, the heating efficiencies should be varied with radii (Yelle 2004). However, the 
Equation (11) of Garcia Mufioz (2007) expresses that they assumed a net heating efficiency 
of 100% for all single photons. Our calculations show also that the ratio of mass loss rates 
decrease to 2.58, roughly equal with the ratio of incident stellar flux, if the net heating 
efficiency is assumed to 100%. This hints that the heating efficiency plays an important role 
in simulating the process of the particle escape. Thus, the realistic heating efficiencies of 
photons in the atmosphere of extra-solar planets should be investigated further by detailed 
radiative transfer calculations. 

An interesting issue is whether a single photon of 20 ev as done in the models of 
Murray-Clay et al. (2009) and Paper I can denote the stellar EUV irradiation. To test 
this assumption, we recalculated the model with a single photon of 20 ev (here we assumed 
that the EUV flux is the same, but we only calculated the radiative transfer of photon of 
20 ev) and found that the mass loss rate is 6.1 x 10^° g which is a factor 2 greater than 
the value predicted by using whole stellar EUV flux. Moreover, our calculations found that 
the maximum mass loss rate (9.55 x 10^° g s~^) is produced when the photon energy is 
about 26 ev, but photons with higher or lower energy will result in lower mass loss rates. 
Murray-Clay et al. (2009) predicted a mass loss rate of 3.3x 10^° g s"-^ when incident flux 
is decreased to F[/y=450 erg s~^ cm~^ integrated on the the solar flux between photon 
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energies of 13.6 ev to 40 ev). Applying the similar stellar flux (F[7y=540 erg s~^ cm~^, the 
value is half of the whole incident flux ), we found that the results calculated by the whole 
EUV flux can be matched well with a single photon of 22ev. By comparing the ionization, 
the location of sonic point, temperature structure and mass loss rates, we found that the 
effect of the whole EUV irradiation can not be substituted by a single photon unless the 
EUV flux is decreased a factor of 2. 

Thus, we can summarize that the effect of the whole EUV irradiation cannot be 
represented simply by a single photon because different photons penetrate to different 
depth and supply different heating. The approximation of a single photon erases the 
characteristics of different photons. In order to obtain reliable results, the contributions 
from all photons should be included fully. 

4. THE TWO-DIMENSIONAL MODELS 

4.1. The effect of tidal force 

In this section we solved the two-dimensional hydrodynamic equations. To separately 
discuss the effect of the tidal force, we assumed the case of non-tidal locking, thus the 
one-dimensional radiative transfer is used (Equation (17)). However, the opacities on the 
day and night sides may be very different. We will further discuss it in Section 4.2. 

The left panels of Figures 3 and 4 show the images of the density (with velocity 
vectors) and temperature. Note that the density is the sum of hydrogen atoms and ions. 
In the 2D model, the number density of hydrogen decreases but the ionization fraction 
increases with radius as does one-dimensional hydrodynamic model. The temperature rises 
with an increase in radius up to r/Rp = 2. The two-dimensional models show that wind 
temperature can attain lOOOOK at about 1.5Rp and then drop to 3000K. 
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The mass loss rates are defined as 

M — 27rm-ii(rjnax)^[ nhUrSin9d9+ / UpUrSinOdO], (22) 

Jit J -k 

where rmax is the upper boundary. We obtained a mass loss rate of M = 1.72 x lO^'^ g 
s~^ which is a factor of 2 smaller than that calculated by one-dimensional hydrodynamics 
model. 

Prom Figures 3 and 4, a prominent characteristic is that the contours of the density and 
temperature are nearly axis-symmetric. It is clear from the left panel of Figure 3 that the 
pattern of flows are non-radial. The flow vectors show curve toward the equator, and the 
radial velocities and densities are very low near the polar regions. We have demonstrated 
this in Figure 5 in which ^"Jf^^^ and ^"''g^Q-) at T=7Rp are showed as a function of ^. It is 
clear that the density and radial velocity decrease rapidly with an increase in 9. The values 
of n^ff^o) ~ ^-^^ u^le=o) ~ ^"^ ^=90° reflect the fact that most particles escape the 
planet from the zones of low latitude. In fact, the zones of middle-high latitude(45°-135'^) 
only contribute to 20% of the mass loss rate due to their low density and radial velocity. 
As a consequence, the mass loss rate of the two-dimensional hydrodynamics is only half of 
that of one-dimensional hydrodynamics. Figure 6 further shows the effect of tidal forces at 
different radii. Due to strong gas pressure, the radial velocities are evidently higher than 
the meridional velocities at the base of the wind (r/Rp=2). However, the ratios of ug/ur 
rapidly increase with an increase in radius. At r/Rp=4, the ratio can attain 0.6 at ^=70°, 
and the value even approaches 1 at r/Rp=7. The trend shows how the tidal forces enforce 
the flows moving in horizonal direction when the gas pressures decrease with an increase in 
radius. In addition, the peaks of ue/ur are almost emerge at ^=70° rather than ^=90° can 
be attributed to the balance of orbital motions and gravity accelerations of host star. 

To further explore the effect of tidal force, we neglected the tidal force and recalculated 
the model. The right panels of Figures 3 and 4 show spherically symmetric distributions 
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for density, velocity and temperature. Comparing the left and right panels of Figures 3 and 
4, we can find that the tidal forces not only supply strong accelerations in direction but 
also decrease the radial velocity in the polar regions. Therefore, we can draw the conclusion 
that the tidal forces not only impel the non-radial flows in the planetary atmosphere but 
also inhibit the escape of particles from the polar regions so that most particles escape the 
bound of the planet from the zones of low latitude. 

4.2. The Case of Tidal Locking 

In Section 4.1, we discussed the case of non-tidal locking. However, about 25% extra- 
solar planets are very close to their host stars thus they should be tidally locked. To inspect 
the case, we recalculated the two-dimensional planetary winds by using two-dimensional 
radiative transfer (Equation (18)), which appropriately depicts the irradiation of the host 
star. 

Figure 7 displays that the decline of density on the night sides with the radius is 
slightly slower than that on day sides. At large radii, the densities of polar regions are 
lower than that of the regions of low-middle latitude. Close to the planet, the meridional 
velocities dominant the velocity field and transfer the mass and energy from the day side 
to the night side (see the left panel of Figure 7 for details). The appearance of strong 
non-radial flows at the base of the wind can be explained by the variations of the optical 
depth. The optical depth of single photon of 20 ev for two-dimensional radiative transfer 
is shown in right panel of Figure 8. The contours of the optical depth extend toward a 
larger radius from the day side to the night side. When the impact factor q approaches 
the radius of the planet, the optical depth is in the order of magnitude of 10"^ — 10^. The 
heating is proportional to e~'^ (Equation (15)), therefore, more energy can be received on 
the day side. With an increase in radius, the radial velocities increase gradually. Finally, 
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the pattern of flows at the base of the wind is nearly in the direction of meridian. In the 
upper atmosphere, the non-radial flows are comparable with the radial flows, but the radial 
velocities are stronger in the zones of low-latitude. 

The distributions of temperature are also non-spherically symmetric due to the 
influence of the optical depth (the left panel of Figure 8). Evidently, the distributions of 
temperature in the night side are controlled by the optical depth. The radiation of the star 
cannot transmit to those regions where the optical depths are greater than unity. Thus 
those regions of lower temperature are located at the proximity of impact factor q~ 1 
(lower left of right panel of Figure 8). The hottest regions is located at 60° < < 180° 
at r/Rp = 1.3 — 1.8, which is caused by the balance of PdV work (PV • v), heating, and 
cooling. To explain the variations of temperature, the contributions of energy at r/Rp=1.35 
by PdV work, heating, cooling and work Wg^t done by external forces (gravity of the 
planet and tidal forces) are shown in Figure 9. Heating from the irradiation and PdV work 
continuously decline to zero from the day side to night side. However, the work done by 
external forces continuously rises with the increase of 6. Due to de = dQ — PdV + Wext, this 
hints that external forces can do positive work due to the contraction of gas, which results 
in the increase of gas temperature (In generally, PdV work expands the atmosphere and 
lowers the gas temperature.). Therefore, the cause of temperature increase from the day to 
night side can be explained by the increase of net heating rate. While 105° ^ 9 < 125°, 
the temperature slightly declines with the decrease of net heating rate. And eventually, the 
temperature becomes a constant because the net heating is maintained at a stable level. 

We also display the angular shces of the density and velocity at the different radii in 
Figures 10-12. For comparison, the winds of one-dimensional radiative transfer are also 
included (solid lines). Seen from Figure 10, the difference in the number density between 
one- and two-dimensional cases is clear. The total number densities of the two-dimensional 
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radiative transfer are much less than those of one-dimensional radiative transfer, by a factor 
of 2-10 (depend on 9). Thus, the mass loss rate also decreases a factor of 4, only 4.3 x 10^ g 
s-\ 

The changes of radial velocity with the angles demonstrate different features between 
the upper and lower atmospheres. In the upper atmosphere (r/Rp = 4 and 7), the radial 
velocities dechne with an increase in angle but rise at ^ = 90° again (the left panel of 
Figure 11). The minimum and maximum radial velocities at r/Rp = 7 are A (9 90°) and 
21 (^ ^ 0°) km s^^, respectively. In comparing with the case of one- dimensional radiative 
transfer, the differences decrease with an increase of radius. 

However, we also note that the radial velocities at r/Rp = 2 almost decline to zero with 
an increase of 9 in the case of two-dimensional radiative transfer. Moreover, the complex 
behaviors at r/Rp =1.5 reflect the difference between two- with one-dimension radiative 
transfer (the left panel of Figure 11). This can be attributed to the effect of optical depth. 
At r/Rp a; 1.5 and 2, the optical depths can vary few order of magnitude from the day side 
to the night side, thus the processes of radiative transfer play an important role. 

In Figure 12, interesting trends in the angular velocity can be found. First, for both 
one- and two- dimensional radiative transfers, the meridional velocities show the sine-like 
profile in the middle-upper atmosphere (the left of Figure 12) and are always negative 
(positive) if 9 is smaller (greater) than 7r/2. In contrast, the meridional velocities are 
always positive at small radii for the case of two-dimensional radiative transfer (the right 
panel of Figure 12), which reflects the transfer of mass and energy from the day side to 
the night side. Second, the meridional velocities are comparable with (or even large than) 
the radial velocities in the middle-high latitude regions of the upper atmosphere. Third, 
the meridional velocities dominant the velocity fleld in the base of the wind for the case of 
two-dimensional radiative transfer, but we have not found similar phenomena in the case of 
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of one-dimensional radiative transfer. 

Therefore, from Figures 7-12, we can conclude that the two-dimensional radiative 
transfer prominently affects the properties of flow. The influence of velocity fleld is 
remarkable at the base of the wind. Due to the large optical depth on the night side, the 
density decrease a factor of few so that the mass loss rate is a factor of 4 lower than that of 
one-dimensional radiative transfer. However, the atmosphere becomes transparent with an 
increase in the radius so that there is moderate difference in the upper atmosphere between 
the one- and two-dimensional radiative transfer. 

5. DISCUSSIONS 

We extended our calculations to different cases. The mass loss rates predicted by our 
models are summarized in Table 1. 

5.1. Boundeiry conditions 

So far, our models used Pr^—I dyn cm to calculate the mass loss from the surface of 
the planet. It is unclear what the suitable lower boundary is? Here we have also computed 
a model in which the lower boundary is set to Pji^=10 dyn cm As expected, the mass 
loss rates are higher than those in the cases of Pr^—I dyn cm Specially, the mass loss 
rates enhance 40-70% for the all cases (see Table 1.). This hints that the mass loss rates 
are very sensitive to lower boundary conditions, and same conclusion has also been found 
by Penz et al. (2008). 

In addition, a possible cause that can lead to the reduction of mass loss rate is the 
day-night temperature contrasts at the lower boundary. Theoretical studies have shown 
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that the day-night temperature contrasts (300K or more) are hkely to occur in the planetary 
photospheres (Showman & Guillot 2002). These predictions are also verified by the observed 
surface temperature distribution of HD 189733b (Knutson et al. 2007). 

In the above calculations, we used the same temperature over the lower boundary. 
In order to further discuss the influence of day-night temperature contrasts, we set the 
temperature at the lower boundary to 



[ 1500i^max(0.5,sine), {9 > 90°). 

Our results did not show more significant anisotropy for the winds (These results are 
not presented in this paper.), and the predicted mass loss rates, 4.05x10^ g s~^, is slightly 
lower than that of model 9 (see Table 1). 



The Swift XRT light curve of HD 189733 shows that the escaping hydrogen is observed 
when the star exhibits a bright fiare that occurrs about 8 hours before the planetary transit 
(Lecavelier des Etangs et al. 2012), which implies a possible correlation between the high 
level of stellar activity and the atmospheric escape of hydrogen. Thus, we also calculated 
the cases with high level of stellar activity (with a solar proxy P10.7 = 200). For the high 
activity, the whole EUV fiux is about 2.5 times of the low activity. Our calculations show 
that the mass loss rates produced by high activity are a factor of 5-10 higher than those of 
low activity, and the temperatures predicted are also higher. By comparing the structures of 
wind, we found that the model with higher stellar activity predicted higher radial velocity, 
and the ratios between meridional and radial velocity are lower. Evidently, higher activities 
supply more UV radiations, especially in the band of higher energy (note that the increases 
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5.2. Dependence of M on Stellar Activities 
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of photon energy are not linear with frequency between the low and high stellar activity). 
More high energy photons mean that the deeper base of the atmosphere can be heated, 
and the hydrodynamic escape parameter, A = {n is mean mass per particle), can 

decrease to a lower value at deeper atmosphere. Thus gas pressure drives a faster expanding 
atmosphere and partly counteracts the effect of the optical depth. It explains the reason 
why the ratio of mass loss rate is only 2.4 for model 8 and 12. 

5.3. Comparison with Observations 

An interesting issue is which mass loss rate predicted by our model is supported by the 
observations. Vidal-Madjar et al. (2003) attributed the excess absorption of HD 209458 
in Ly q; to a mass loss of 10^° g s~^. Subsequent ID models also support this. However, 
the influences of two-dimension and charge exchange between the stellar wind and the 
planetary escaping exosphere were neglected by all the ID models. Recently, Lecavelier des 
Etangs et al. (2012) updated the mass loss rate to 10^ g s~^ for the atmospheric escape of 
HD 189733b though the previous mass loss rate requested in their model is in the order 
of magnitude of 10^- 10^^ g s~^. Moreover, Holmstrom et al. (2008) found that a mass 
loss rate of 7 x 10^ g s"^ can explain the observations of HD 209458b when the charge 
exchange is included. We must be careful in the interpretation of mass loss rates because 
only atomic hydrogen contributes to the excess absorption of Ly a. In our models, the 
mass loss rates of atomic hydrogen are 3.78 x 10^ g s~^, 9.85 x 10^ g s~^, 7.60 x 10^ g s~^ 
and 1.45 x 10^° g s~^ for models 9-12, respectively. If the charge exchange indeed results 
in the excess absorption of Ly a, our results favor the low activity for HD 209458 when it 
is transited by HD 209458b. However, the possibility that planetary origin hydrogen result 
in excess absorption of Ly a cannot be excluded because Koskinen et al. (2010) explained 
this by solely absorption in the upper atmosphere based on the ID assumption. Their 
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model denotes the averaged property of the atmosphere, but is very different with our 2D 
hydrodynamic model with consistent radiative transfer. The differences appear in several 
aspects. First, they estimated a mass loss rate of 1-10x10^" g s~^, which is higher than our 
result. Second, in order to fit the observations, they estimated a number density of 2.6 x 
10^ cm~^ for hydrogen atoms at r=2.9Rp. However, we found that the number densities 
of hydrogen atoms in our model are about 10^-10^ cm~'^ at r=2.9Rp (depend on 6), and 
the number density distributions of hydrogen atoms show clear variations from the day-side 
to night-side (see Figure. 13). Third, they assumed an ionization of 100% at r=2.9Rp. 
Our results show that the ionization structures of hydrogen strongly depend on 9 for the 
case of tidal-locking. For example, ^=0.16, 0.19 and 1.08 at 9=0, 7t/2 and tt at r=2.9Rp. 
However, our models predicted the similar temperature distributions (Tr^TOOO-IIOOOK) 
below r=2.9Rp as estimated by their models. Thus, one-dimensional models are not 
accurate enough in calculating the HI transit depth. Our results can be applied further to 
check if the upper atmosphere contributes significant absorption. 

6. CONCLUSIONS 

HST observations of HD 209458b and HD 189733b in the UV band have revealed 
the signatures of atmosphere escape. Theoretical progress is following the observations, 
for example, reasonable mass loss rates have been predicted by several one-dimensional 
hydrodynamic models (Yelle 2004; Tian et al. 2005; Garcia Munoz 2007; Penz 2008; 
Murray-Clay et al. 2009; Guo 2011). Two-dimensional hydrodynamic model have also 
revealed a possible anisotropy in the planetary wind (Stone & Proge 2009). However, there 
exists still a few key problems to be solved, such as the proper descriptions to radiative 
transfer and charge-exchange reactions between the planetary and stellar wind particles. 
They play an important role in fully explaining the observations. 
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Table 1: 


The mass loss rates for HD 209458b. 
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In this paper, we presented a two-dimensional multi-fluid hydrodynamic model by 
incorporating the two-dimensional radiative transfer. Prom this work, our conclusions are 
the following. 

1. In the assumption of non-tidal locking, the planetary winds are not spherically 
symmetric. Tidal forces can result in significant horizontal movements in the planetary 
winds and lower density and radial velocity near the polar regions so that the mass is 
almost lost fully from the low-latitude zones. 

2. In the case of tidal locking, we found that the wind is controlled by meridional 
velocity in the regions of r/Rp ~ 1.5. The large meridional velocities transfer the mass and 
energy from the day side to the night side at the base of the wind. But the extent decreases 
with an increase in the radius due to the decrease of the optical depth. At large radii, the 
wind is siginificant difference with that of one-dimensional radiative transfer. Due to the 
dechnes of density, the mass loss rate is also lower than that of non-tidal locking. 

3. The mass loss rates are affected by both two-dimensional radiative transfer and 
hydrodynamics. Comparing with the results of one-dimensional hydrodynamic models, the 
models of two-dimensional hydrodynamics predicted the mass loss rates of 1.72 and 0.43 x 
10^° g s~^ for one-dimensional and two-dimensional radiative transfer, which are a factor of 
2 and 7 smaller than that of ID hydrodynamic model. 
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Fig. 1. — Schematic diagram of the calculation method to the optical depth. The circles 
represent the nodal points at other radii, i represents the radial grid, j denotes the angular 
grid. 




Fig. 2. — The one-dimensional wind model for HD 209458b. Number densities (upper left), 
velocities (upper right), temperatures (lower left) and the ionization fraction (lower right) 
are plotted as the functions of altitude. 
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Fig. 3. — The results of two-dimensional hydro dynamic model with one-dimensional radiative 
transfer. Left: density and velocity vectors for the case with tidal forces. Right: density and 
velocity vectors for the case without tidal forces. The host star is located toward the right 
(corresponding to 6^ = 0). 
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Fig. 4. — The results of two-dimensional hydrodynamic model with one-dimensional radiative 
transfer. Left: temperature distributions for the case with tidal forces. Right: temperature 
distributions for the case without tidal forces. The host star is located toward the right 
(corresponding to 6 = 0). 
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Fig. 5. — Left: the radial velocity as a function of 9 at r/Rp=7. Right: the distributions 
of number density as a function of 9 at r/Rj,=7. 




Fig. 6. — The values of ug/ur as a function of 9 at r/Rp=2, 4 and 7. 
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Fig. 7. — The results of the two-dimensional hydro dynamic model with two-dimensional 
radiative transfer. Left: density and velocity vectors at 1 ^ r/Rp ^ 7. Right: density and 
velocity vectors at the atmospheric base (1 ^ r/Rp ^ 2.6). 
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Fig. 8. — The results of the two-dimensional hydrodynamic model with two-dimensional 
radiative transfer. Left: the distributions of temperature. Right: the distributions of the 
optical depth. In the right panel the optical depth lager than 10 is marked with "dark". 
The "gray" region denotes the case of 10^^ < r < 1. In the transparent region, the optical 
depth is smaller than 10~^. 
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Fig. 9. — Temperatures and heating rates (unit mass) in the model. In the lower panel, 
the solid line represents the net heating rate in the atmosphere. The dashed, dotted and 
dash-dotted hues denote the PdV work, heating from host star and coohng. The term MVext 
(dash-dot-dotted) represents the work made by the gravity of planet and tidal forces. 
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Fig. 10. — Number densities at r/Rp =1.5, 2, 4 and 7 (from top to bottom). Solid lines 
represent the case of one-dimensional radiative transfer. Dashed lines represent the case of 
two-dimensional radiative transfer. Note: The number densities are the sum of hydrogen 
atoms and ions. 
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Fig. 11. — Left: radial velocities at r/Rp =2, 4 and 7. Right: radial velocities at r/Rp =1.5. 
Solid lines represent the case of one-dimensional radiative transfer. Dashed lines represent 
the case of two-dimensional radiative transfer. 
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Fig. 12. — Left: meridional velocities at r/Rp =2, 4 and 7. Right: meridional velocities at 
r/Rp =1.5. Solid lines represent the case of one-dimensional radiative transfer. Dashed lines 
represent the case of two-dimensional radiative transfer. 
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Fig. 13. — The number density distributions of hydrogen atoms in the case of the two- 
dimensional hydrodynamic model with two-dimensional radiative transfer. 
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